***************
* This script produces tables and figures based on the estimates of the part-time penalty by gender and isco 2-digit occupation
* Author: Daniel Kopp
***************

	clear
	clear matrix
	clear mata
	set more off	
	set maxvar 8000

	use "Misc_files\share_seeking_fulltime_gender_isco2.dta"

	cap gen y = 1
	cap gen x = 1	
	reg y x

	nlcom (base0: _b[x]) (base1: _b[x]), post
	estimates store base0
	
	******************************************************************	
	******************************************************************
	* Part-time penalty (for men and women separately) by isco 2
	******************************************************************
	******************************************************************

	foreach x in inter_ratio {
		cap gen penalty_`x' = . 
		cap gen t_penalty_`x' =. 
		cap gen se_penalty_`x' =.
		cap gen lo_penalty_`x' = .
		cap gen hi_penalty_`x' = .				
	}

	eststo clear
	estimates use "$save_path\workvolume_by_isco2_int"
	estimates store part2_i2_interact	
	
	estimates restore part2_i2_interact
		local mean_sample = e(mean_sample)
		local mean_disp = e(mean_disp)		
		local N = e(N)
		matrix b = e(b) 
		matrix b = b[1,1..120]
		matrix V = e(V)
		matrix V = V[1..120,1..120]
		local obs : disp %9.0f `N'
		quietly ereturn post b V, dep(contact_button_clicked) obs(`obs') 	
			estadd local mean_sample   " `mean_sample' "
			estadd local mean_disp   " `mean_disp' "			
			estadd local obs  " `obs' "															
	eststo part2_i2_interact_sm	
		
	esttab part2_i2_interact_sm	, keep(isco2_11Xparttime isco2_12Xparttime isco2_13Xparttime isco2_14Xparttime isco2_21Xparttime isco2_22Xparttime isco2_23Xparttime isco2_24Xparttime isco2_25Xparttime isco2_26Xparttime isco2_31Xparttime isco2_32Xparttime isco2_33Xparttime isco2_34Xparttime isco2_35Xparttime isco2_41Xparttime isco2_42Xparttime isco2_43Xparttime isco2_44Xparttime isco2_51Xparttime isco2_52Xparttime isco2_53Xparttime isco2_54Xparttime isco2_61Xparttime isco2_62Xparttime isco2_71Xparttime isco2_72Xparttime isco2_73Xparttime isco2_74Xparttime isco2_75Xparttime isco2_81Xparttime isco2_82Xparttime isco2_83Xparttime isco2_91Xparttime isco2_92Xparttime isco2_93Xparttime isco2_95Xparttime isco2_96Xparttime isco2_99Xparttime isco2_999Xparttime)
	* One coefficient is (insignificantly) positive (isco2_22)
				
	* Gen number of observations and mean contact rate by isco-2 occupation
	gen no_obs_isco = . 
	gen mean_contact_isco = . 
	foreach i in 11 12 13 14 21 22 23 24 25 26 31 32 33 34 35 41 42 43 44 51 52 53 54 61 62 71 72 73 74 75 81 82 83  91 92 93 96 99  {		
		estimates use "$save_path\summary_isco`i'"
		matrix list e(mean)
		matrix c = e(mean)
		local obs_`i' = e(N)
		global mean_`i' =  c[1,1]
		disp "`obs_`i''"
		disp "${mean_`i'}"
		replace no_obs_isco = `obs_`i'' if s_isco_08_2==`i' 
		replace mean_contact_isco = ${mean_`i'} if s_isco_08_2==`i' 
	}	
		
	* Without  95 since there the parttime coefficient could not be estimated
	foreach i in 11 12 13 14 21 22 23 24 25 26 31 32 33 34 35 41 42 43 44 51 52 53 54 61 62 71 72 73 74 75 81 82 83  91 92 93 96 99  {		
	
	* Scale coefficients
	estimates restore part2_i2_interact_sm		
			local scale_`i': disp %9.6g  100/${mean_`i'}
			local N = e(N)
			matrix b = e(b) * `scale_`i''
			matrix V = e(V) * `scale_`i'' ^2
			local obs : disp %8.0f `N'
			quietly ereturn post b V, dep(contact_button_clicked) obs(`obs') 								
	eststo 	part2_i2_interact_sc	
			
	* Ratio of female to male part-time penalty 
	estimates restore part2_i2_interact_sc
	local mean_disp = e(mean_disp)		
	nlcom		(parttimeXgeschlecht: (_b[isco2_`i'Xparttime] + _b[isco2_`i'XparttimeXgeschlecht])/_b[isco2_`i'Xparttime]) , post
	estadd local mean_disp   " `mean_disp' "				
	estimates 	store inter_ratio_i2_`i'	
			
	}
	
	* Gen variable with ratio
	* Without  95 since there the parttime coefficient could not be estimated
	foreach i in 11 12 13 14 21 22 23 24 25 26 31 32 33 34 35 41 42 43 44 51 52 53 54 61 62 71 72 73 74 75 81 82 83  91 92 93 96 99  {		
		foreach x in  inter_ratio  {
		
		estimates restore `x'_i2_`i'	
			cap	local 	penalty_`x' 		= _b[parttimeXgeschlecht] 
			cap	replace penalty_`x' 	 	= `penalty_`x'' 		if s_isco_08_2==`i'		
			cap local 	t_penalty_`x' 	 	= abs(_b[parttimeXgeschlecht] / _se[parttimeXgeschlecht] )
			cap replace t_penalty_`x' 	 	= `t_penalty_`x'' 		if s_isco_08_2==`i'		
			cap local 	se_penalty_`x' 	 	=  _se[parttimeXgeschlecht] 
			cap replace se_penalty_`x' 	 	= `se_penalty_`x'' 		if s_isco_08_2==`i'		
			cap local 	lo_penalty_`x' 		=  _b[parttimeXgeschlecht] - 1.96 * _se[parttimeXgeschlecht] 
			cap replace lo_penalty_`x' 		= `lo_penalty_`x'' 		if s_isco_08_2==`i'		
			cap local 	hi_penalty_`x' 		=  _b[parttimeXgeschlecht] + 1.96 * _se[parttimeXgeschlecht] 
			cap replace hi_penalty_`x' 		= `hi_penalty_`x'' 		if s_isco_08_2==`i'			
		estimates drop  `x'_i2_`i'	
		}
	}
	
	
	******************************
	* Heterogeneity of the ratio of the female to male part-time penalty by ISCO 2-digit occupations (with and without shrinkage)
	******************************
		
	* The following two "ratios" are difficult to interpret: 
		* Isco 22 (Health professionals): Male part-time coefficient positive and female one negative
			* -> the resulting ratio is negative
		* Isco 11 (Chief executives): Male part-time coefficient very negative and female one positive 
			* -> the resulting ratio is negative	
	
	foreach coef in penalty_inter_ratio  {
	
	****
	* First, the shrunken coefficients:
	****
	
	* For the empirical shrinkage procedure, see "Koedel et al. 2015: Value-Added Modeling: A Review" for the formula
	
	cap drop lambda_j
	gen   lambda_j = se_`coef'^2

	quietly sum lambda_j
	local mean_lambda_j = r(mean)
	
	quietly sum `coef'
	local mean_beta = r(mean)		
	disp "The mean of the unshrunken coefficients is `mean_beta'"
	local sigma_sq = r(Var)			
	disp "`sigma_sq'"								
														 
	cap drop alpha
	gen alpha = `sigma_sq' / (`sigma_sq' + lambda_j)
	
	* Gen EB estimate:
	cap drop `coef'_adj
	gen `coef'_adj = alpha*`coef'+(1-alpha)*`mean_beta'

	sum `coef' `coef'_adj,d

	* Gen se of EB estimate:
	cap drop se_`coef'_adj
	gen se_`coef'_adj = (alpha*lambda_j)^0.5   // See the formula in Herrmann et al. 2016: Shrinkage of Value-Added Estimates, page 3 

	sum  se_`coef' se_`coef'_adj
	
	cap drop t_`coef'_adj 	
	gen t_`coef'_adj 	 	= abs(`coef'_adj / se_`coef'_adj )

	cap drop hi_`coef'_adj lo_`coef'_adj
	gen hi_`coef'_adj  = `coef'_adj  + 1.96 * se_`coef'_adj
	gen lo_`coef'_adj  = `coef'_adj  - 1.96 * se_`coef'_adj
		
	preserve 	
	cap drop not_sign_adj
	gen not_sign_adj = lo_`coef'_adj<1 & hi_`coef'_adj>1 & hi_`coef'_adj<.
	tab not_sign_adj,m

	sort `coef'_adj
	*cap drop n_adj
	drop if s_isco_08_2==22		// see remark above
	drop if s_isco_08_2==11		// see remark above	
	drop if s_isco_08_2==99		// no information
	drop if t_`coef'==.
	
	gen n_adj=_n	
	tab not_sign_adj,m
		
	count if `coef'_adj>1		//  9
	count if `coef'_adj<1		// 27
	* Save the mean
		quietly sum `coef'_adj , d	
		local mean_ratio = r(mean)
	
	labmask n_adj , values(s_isco_08_2) decode
	estpost tabstat `coef'_adj , by(s_isco_08_2) statistics(mean) notot
	eststo shrinkage_ratio_co
	estpost tabstat se_`coef'_adj, by(s_isco_08_2) statistics(mean) notot
	eststo shrinkage_ratio_se
	labmask n_adj , values(s_isco_08_2) 
		
	replace lo_`coef' = -1 if lo_`coef'<-1  // we cap the CI at -1 and +3
	replace hi_`coef' = 3 if hi_`coef'>3  // we cap the CI at -1 and +3		
	gen n_adjoff_m = n_adj - 0.2  

	twoway 	(rcap 		lo_`coef' hi_`coef' 		n_adj 		, lcolor(gs10) msize(vtiny)) || ///
			(rcap 		lo_`coef'_adj hi_`coef'_adj n_adjoff_m 		, lcolor(black*0.85) msize(vtiny)) || ///
			(scatter 	`coef' 						n_adj		, mcolor(gs10) msize(small) )  ||  ///			
			(scatter 	`coef'_adj 					n_adjoff_m		, mcolor(black*0.85) msize(small))  ,  ///
			graphregion(color(white))    ytitle(Ratio of female to male part-time penalty)  ///
			xtitle(ISCO-2 occupations)    xlabel(1(1)36, labsize(small) alt valuelabel) ///
			yline(1, lcolor(gs10) lpattern(dash)) yline(`mean_ratio', lcolor(black) lpattern(dash)) ///			
			legend(order(3 "Original coefficient" 4 "Shrunken coefficient"))	
	graph export "$results_part_time/figure_d8.eps", as(eps) replace
					   
	* Gen global with order of isco occupations for table
	decode n_adj, gen(n_adj_string)
	destring n_adj_string, gen(n_adj_new)
	local order_isco_shrunken ""
	forvalues i = 1/35 {
		local val_lab = n_adj_new[`i']
		local order_isco_shrunken "`order_isco_shrunken' `val_lab'"
	}
	global order_isco_shrunken "`order_isco_shrunken'"
	disp "$order_isco_shrunken"
	
	restore
	
	****
	* Second, the unshrunken coefficients:
	****
	
	preserve 	
	cap drop not_sign
	gen not_sign = lo_`coef'<1 & hi_`coef'>1 & hi_`coef'<.
	tab not_sign,m

	sort `coef'
	*cap drop n
	drop if s_isco_08_2==22		// see remark above
	drop if s_isco_08_2==11		// see remark above	
	drop if s_isco_08_2==99		// no information
	drop if t_`coef'==.
	
	gen n=_n	
	tab not_sign,m
	
	count if `coef'>1		//  
	count if `coef'<1		// 
	* Save the mean
		quietly sum `coef'  , d	
		local mean_ratio = r(mean)
		
	estpost tabstat `coef' , by(s_isco_08_2) statistics(mean) notot
	eststo ratio_co
	estpost tabstat se_`coef', by(s_isco_08_2) statistics(mean) notot
	eststo ratio_se
	restore
	}
	 			
	* Show the shrunken and unshrunken estimates in table
	esttab shrinkage_ratio_co shrinkage_ratio_se ratio_co ratio_se ,  ///
			cells("mean(fmt(%9.2gc))") label nonumbers  varlabels(`e(labels)') varwidth(45)  ///
			order($order_isco_shrunken) ///
			mgroups("\textbf{Shrunken}" "\textbf{Not shrunken}", pattern(1 0 1 0) span prefix(\multicolumn{@span}{c}{) suffix(})) ///
			mtitle("Ratio" "S.E." "Ratio" "S.E." )  collabels(none) ///
			addnote("Source: Seco Jobroom data")				
	esttab shrinkage_ratio_co shrinkage_ratio_se ratio_co ratio_se using "$results_part_time/figure_d8_table.tex", replace ///
			cells("mean(fmt(%9.2gc))") label nonumbers  varlabels(`e(labels)') varwidth(45) frag ///
			order($order_isco_shrunken) ///
			mgroups("\textbf{Shrunken}" "\textbf{Not shrunken}", pattern(1 0 1 0) span prefix(\multicolumn{@span}{c}{) suffix(})) ///
			mtitle("Ratio" "S.E." "Ratio" "S.E." )  collabels(none) ///
			addnote("Source: Seco Jobroom data")					
			

				

	
	
	
